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Abstract 

A class of distributed systems with a cyclic interconnection structure is considered. These systems 
arise in several biochemical applications and they can undergo diffusion driven instability which leads 
to a formation of spatially heterogeneous patterns. In this paper, a class of cyclic systems in which 
addition of diffusion does not have a destabilizing effect is identified. For these systems global stability 
results hold if the "secant" criterion is satisfied. In the linear case, it is shown that the secant condition 
is necessary and sufficient for the existence of a decoupled quadratic Lyapunov function, which extends 
a recent diagonal stability result to partial differential equations. For reaction-diffusion equations with 
nondecrcasing coupling nonlinearities global asymptotic stability of the origin is established. All of the 
derived results remain true for both linear and nonlinear positive diffusion terms. Similar results are 
shown for compartmental systems. 

1 Introduction 

The first gene regulation system to be studied in detail was the one responsible for the control of lactose 
metabolism in E. Coli, the lac operon studied in the classical work of Jacob and Monod [H[2]. Jacob and 
Monod's work led Goodwin [3] and later many others [l[^[6l[7l[8l[^ [M[TTl[T^[T^[Tl[T5] to the mathematical 
study of systems made up of cyclically interconnected genes and gene products. In addition to gene 
regulation networks, cyclic feedback structures have been used as models of certain metabolic pathways |16j . 
of tissue growth regulation [17] , of cellular signaling pathways [18] , and of neuron models [19] . 

Generally, cyclic feedback systems (of arbitrary order) were shown by Mallet-Paret and Smith [20]I21] 
to have behaviors no more complicated that those of second-order systems: for precompact trajectories, 
w-limit sets can only consist of equilibria, limit cycles, or heteroclinic or homoclinic connections, just as 
in the planar Poincare-Bendixson Theorem. When the net effect around the loop is positive, no (stable) 
oscillations are possible, because the overall system is monotone [22]. On the other hand, inhibitory or 
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"negative feedback" loops give rise to the possibility of periodic orbits, and it is then of interest to provide 
conditions for oscillations or lack thereof. 

Besides the scientific and mathematical interest of the study of cyclic negative feedback systems, there 
is an engineering motivation as well, which rises from the field of synthetic biology. Oscillators will be 
fundamental parts of engineered gene bacterial networks, used to provide timing and periodic signals to 
other components. A major experimental effort, pioneered by the construction of the "repressilator" by 
Elowitz and Leibler [23J, is now under way to build reliable oscillators with gene products. Indeed, the 
theory of cyclic feedback systems has been proposed as a way to analyze the repressilator and similar 
systems j2Hl25] . 

In order to evaluate stability properties of negative feedback cyclic systems, [9] and |15| analyzed the 
Jacobian linearization at the equilibrium, which is of the form 

~ -ax • • • -b r 
bx — 02 '•• 
A = 6 2 -03 *'■ : 

: '■■ '•• '■• 

••• 6 n _i -a, 

di > 0, bi > 0, i = 1, • • • ,n, and showed that A is Hurwitz if the following sufficient condition holds: 

61 ' - bn < sec(7r/n)" (2) 
ax ■ ■ ■ a n 

This "secant criterion" is also necessary for stability when the a,'s are identical. 

An application of the secant condition in a "systems biology" context was in Kholodenko's [18] (see 
also |26j) analysis of a simplified model of negative feedback around MAPK (mitogen activated protein 
kinase) cascades. MAPK cascades constitute a highly conserved eukaryotic pathway, responsible for some 
of the most fundamental processes of life such as cell proliferation and growth [2711281(29] . Kholodenko used 
the secant condition to establish conditions for asymptotic stability. 



(1) 



1.1 Global stability considerations 

It appears not to be generally appreciated that (local) stability of the equilibrium in a cyclic negative 
feedback system does not rule out the possibility of periodic orbits. Indeed, the Poincare-Bendixson 
Theorem of Mallet-Paret and Smith |2CHI21| allows such periodic orbits to coexist with stable equilibria. 
As an illustration consider the system 

Xi = -Xx + <p(X3) 

X2 = -X2 + X1 (3) 

X3 = ~X3 + X2 

where 

<p(Xs) = e" 10 ^ 3 " 1 ) + 0.1sat(25( X3 - 1)), (4) 
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and sat(-) := sgn(-) min{l, | • |} is a saturation^] function. The function @ is decreasing, and its slope has 
magnitude 63 = 7.5 at the equilibrium xi = X2 = X3 = 1- With a\ = 02 = 03 = b\ = 62 = 1 and n = 3, the 
secant criterion ([2]) is satisfied and, thus, the equilibrium is asymptotically stable. However, simulations 
in Fig. [1] show the existence of a periodic orbit in addition to this stable equilibrium. 




Figure 1: Trajectory of ([3]) starting from initial condition x = [1-2 1-2 1.2] T , projected onto the xi~X2 
plane. 

To delineate global stability properties of cyclic systems with negative feedback, [30] studied (by building 
on a passivity interpretation of the secant criterion in [31]) the nonlinear model 

xi = -fi(xi) - g n {x n ) 
X% = ~f2{x 2 ) + gi(xi) 

(5) 

Xn — fn\Xn) 9n—l\Xn—l) 

and proved global asymptotic stability of the original under the conditions 

afiia) > 0, a 9i (a) > 0, Vo" G M\{0}, (CI) 



9i{<r) 



< 7i , Va G M\{0}, 



71 • • - 7n < sec(7r/n) 



(C2) 
(C3) 



lim 

|aci| — * 00 



9i{a)da = 00. 



(C4) 



The conditions (CI) (C4) encompass the linear system (d])-© in which fi(xi) = aiXi, gi{xi) = 6jXj, and 
7i = fej/aj. 

A crucial ingredient in the global asymptotic stability proof of [30] is the observation that the secant 
criterion ([2]) is necessary and sufficient for diagonal stability of ([T]), that is for the existence of a diagonal 
matrix D > such that 

A T D + DA < 0. (6) 



1 One can easily modify this example to make ip(-) smooth while retaining the same stability properties. 
2 In the rest of the paper we assume that an equilibrium exists and is unique (see [30] for conditions that guarantee this) 
and that this equilibrium has been shifted to the origin with a change of variables. 
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Using this diagonal stability property, [30] constructs a Lyapunov function for ([5]) which consists of a 
weighted sum of decoupled functions of the form Vi(xi) = f^ 1 gi{&) da. In the linear case this construction 
coincides with the quadratic Lyapunov function V = x T Dx. 

1.2 Spatial localization 

Ordinary differential equation models such as described above implicitly assume that reactions proceed in 
a "well-mixed" environment. However, in cells, certain processes are localized to membranes (activation of 
pathways by receptors), to the nucleus (transcription factor binding to DNA, production of mRNA), to the 
cytoplasm (much of signaling), or to one of the specialized organelles in eukaryotes. The exchange of chem- 
ical species between these spatial domains has been found to be responsible for dynamical behavior, such 
as emergence of oscillations, in fundamental cell signaling pathways, see for instance [32]. These exchanges 
often happen by random movement (diffusion), although transport mechanisms and gated channels are 
sometimes involved as well. 

When each of a finite set of spatial domains is reasonably "well-mixed," so that the concentrations of 
relevant chemicals in each domain are appropriately described by ordinary differential equations (ODEs), 
a compartmental model may be used. In a compartmental model, several copies of an ODE system are 
interconnected by "pipes" that tend to balance species concentrations among connected compartments. 
The overall system is still described by a system of ODEs, but new dynamical properties may emerge from 
this interconnection. For example, two copies of an oscillating system may synchronize, or two multi-stable 
systems may converge to the same steady state. 

On the other hand, if a well-mixed assumption in each of a finite number of compartments is not 
reasonable, a more appropriate mathematical formalism is that of reaction-diffusion partial differential 
equations (PDEs) [33. 34, 35, 36, 37]: instead of a dynamics x = f(x), one considers equations of the general 
form 

| = DA, + /(,), g=0, (7) 

where now the vector x = x(£, t) depends on both time t and space variables £ belonging to some domain 
ft, Ax is the Laplacian of the vector x with respect to the space variables, D is a matrix of positive diffusion 
constants, and dx/dv denotes the directional derivative in the direction of the normal to the boundary 
dVt of the domain fi, representing a no-flux or Neumann boundary condition. (Technical details are given 
later, including generalizations to more general elliptic operators that model space-dependent diffusions.) 

Diffusion plays a role in generating new behaviors for the PDE as compared to the original ODE 
x = f{x). In fact, one of the main areas of research in mathematical biology concerns the phenomenon of 
diffusive instability, which constitutes the basis of Turing's mechanism for pattern formation [38 , 39 , 40ll41j , 
and which amounts to the emergence of stable non-homogeneous in space solutions of a reaction-diffusion 
PDE. The Turing phenomenon has a simple analog, and is easiest to understand intuitively, for an ODE 
consisting of two identical compartments |41}I42|. Also in the context of cell signaling, and in particular 
for the MAPK pathway mentioned earlier, reaction-diffusion PDE models play an important role [43j . 

If diffusion coefficients are very large, diffusion effects may be ignored in modeling. As an illustration, 
the stability of uniform steady states is unchanged provided that the diffusion coefficient D is sufficiently 
large compared to the "steepness" of the reaction term /, measured for instance by an upper bound a on 
its Lipschitz constant or equivalently the maximum of its Jacobians at all points (for chemical reaction 
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networks, this is interpreted as the inverse of the kinetic relaxation time, for steady states). Introducing an 
energy function using the integral of \dx/d^\ 2 , and then integrating by parts and using Poincare's inequality, 
one obtains an exponential decrease of this energy, controlled by the difference of a and D ( [UJ , Chapter 
11). For instance, Othmer [35] provides a condition Dfi > a in terms of the smallest nonzero eigenvalue of 
the Neumann Laplacian {Ax + fix = 0, £ £ £1; dx/dv = 0, £ £ 957} to guarantee exponential convergence 
to zero of spatial nonuniformities, and estimates that his condition is met for intervals f2 = [0, L] of length 
L ~ lOfim, with diffusion of at least about 4x 10 -8 cm 2 /sec and a ~ 10 _1 sec. 

On the other hand, if diffusion is not dominant, it is necessary to explicitly incorporate spatial inho- 
mogeneity, whether through compartmental or PDE models. The goal of this paper is to extend the linear 
and nonlinear secant condition to such compartmental and PDE models, using a passivity-based approach. 
To illustrate why spatial behavior may lead to interesting new phenomena even for cyclic negative feedback 
systems, we take a two-compartment version of the system shown in J3]): 

Xi = -Xi + f(X3) + D (Vi ~ Xi) 

X2 = -X2 + xi + D(m - X2) 

X3 = ~X3 +X2 + D{r}3 - X3) 

m = -m + ^>{m) + D (xi -m) (8) 

m = -m + m + D (x2 - m) 

m = -m + m + D (x3 - m) 

and pick D = 10~ 4 . We simulated this system with initial condition [1.4945 1.3844 1.0877 1 1 1] T , so 
that the first-compartment Xi(0) coordinates start approximately on the limit cycle, and the second- 
compartment r]i(0) coordinates start at the equilibrium. The resulting simulation shows that a new os- 
cillation appears, in which both components oscillate, out of phase (no synchronization), with roughly 
equal period but very different amplitudes. Figure [2] shows the solution coordinates xi an d Vi plotted 

1.6 
1.4 
1.2 



0.8 



"1982 1984 1986 1988 1990 1992 1994 1996 1998 2000 
t 

Figure 2: New oscillations in two-compartment system: xi (solid) and r/i (dashed) shown. 

on a window after a transient behavior. This oscillation is an emergent behavior of the compartmental 
system, and is different from the limit cycle in the original three-dimensional system. (One may analyze 
the existence and stability of these orbits using an ISS-like small-gain theorem.) 
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Our goal is to show that -in contrast to this example-if the secant condition does apply to a negative 
cyclic feedback system, then no non-homogeneous limit behavior can arise, in compartmental or in PDE 
models, no matter what is the magnitude of the diffusion effect. 



2 Problem formulation 



In this paper we extend the linear and nonlinear results of [9, 15 , 30J to spatially distributed models that 
consist of a cyclic interconnection of n reaction-diffusion equations 

fat = V-(hi(fa)Vfa) - h(fa) - g n (fa) 
fat = V- (h 2 (fa) Vfa) - f 2 (fa) + 9i(fa) 

. (RD) 

fat = V • (h n (lp n )Vfa) - f n (fa) + 9n-l(fa-l) 

where fa denotes the state of the ith subsystem which depends on spatial coordinate £ and time t, fa(£, t), 
and fi, gi, hi denote static nonlinear functions of their arguments. We consider a situation in which the 
spatial coordinate £ := (£i, . . . , £ r ) belongs to a bounded domain f2 in W, r = 1, 2 or 3, with a smooth 
boundary d£l and outward unit normal v. The state of each subsystem satisfies the Neumann boundary 
conditions, dfa/du := fa v = on dQ, V^j is the gradient of ipi, V • v is the divergence of a vector v, and 
the domain of the r-dimensional Laplacian A := V • V is given by [461 147] 

V(A) := {ipi G H 2 (Q), ip iU = on dfl} . (DM) 

Here, H 2 (Q) denotes a Sobolev space of square integrable functions with square integrable second distri- 
butional derivatives. The standard L^iP) inner product is given by 

:= / ^ T (OttO<% 
Jn 

rp 

where d£ := d^i • • • d£ r and ip := \ fa ■ ■ ■ tp n ] . 

As explained in the introduction, the study of stability properties for distributed system |(RD)| is 
important in many biological applications. Our first result, presented in Section [31 studies the linearization 
of |(RD)| and shows that the secant condition ^ is sufficient for the exponential stability despite the 
presence of diffusion terms. It further shows that the secant condition is necessary and sufficient for the 
existence of a decoupled Lyapunov function, thus extending the diagonal stability result of [30] to partial 
differential equations. The next result of the paper, presented in Section [H studies the nonlinear reaction- 



diffusion equation (RD) and proves global asymptotic stability of ip = under assumptions that mimic 
the conditions (C1)|(C3) of [30], and under the additional assumptions that the functions gi(-) and hi(-), 
i = 1, • • • ,n, be nondecreasing and positive, respectively. This additional assumption on the ^-functions 
ensures convexity of the Lyapunov function which is a crucial property for our stability proof. Indeed, 
a similar convexity assumption has been employed in [48] to preserve stability in the presence of linear 
diffusion terms. Finally, Section [5] studies a compartmental ordinary differential equation model instead of 
the partial differential equation |(RD)| and proves global asymptotic stability using the same nondecreasing 
assumption for g^s. 
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3 Cyclic interconnection of linear reaction-diffusion equations 



We start our analysis by considering an interconnection of spatially distributed systems |(RD)| with 

fi(i/ji) := aitpi, gi(tpi) := htpt, hi(ipi) := a, i = 1, . . . , n (9) 



where each m, b{, and c% represents a positive parameter. In this case, system ["(RD)| simplifies to a cascade 
connection of linear reaction-diffusion equations where the output of the last subsystem is brought to the 
input of the first subsystem through a negative unity feedback. Abstractly, the dynamics of system |(RD)| - 
(DM) with fi(-), gi(-), and /ij(-) satisfying Q are given by 



^ = Aip := CAip + A ip 



where denotes the vector Laplacian, that is Aijj := [ Ai/j\ 
0, and 



AV>n ] T , C : = diag{[ C\ 





-Oi 





-b n 










h 


- a 2 
















b 2 




, Oi > 0, I 


H > 0, i = 1, . . 


. ,n 
























b n -i - a n 









(LRD) 

c„ ]} > 



3.1 Exponential stability and the secant criterion in one spatial dimension 

In this section, we focus on systems with one spatial dimension £ G £1 := (0,1). We show that operator 



A with (DM) generates an exponentially stable strongly continuous (C Q ) semigroup T(t) on L^iO, 1) if 
the secant criterion ([2]) is satisfied. We note that the exponential stability of T(t) in Theorem [1] can be 
also established using a Lyapunov based approach that we develop for systems with two or three spatial 
coordinates. However, the proof of Theorem[T]is of independent interest because of the explicit construction 



of the C -semigroup and block-diagonalization of operator (LRD) (DM) (which is well suited for a modal 
interpretation of stability results in one spatial coordinate) . 

It is well known (see, for example |47| ) that the operator with Neumann boundary conditions is 
self-adjoint with the following set of eigenfunctions {<fik} and corresponding eigenvalues {f/t}: 

MO = 1, ¥>i(£) = v^cos l7r£, I £ N, 
z/ = 0, vi = -(In) 2 , I e N. 

Since the eigenfunctions {<Pk} represent an orthonormal basis of £2(0, 1) each ipi(£,,t) can be represented 
as 



fc = 



where Xi k(t) denote the spectral coefficients given by 



Xi,k(t) = {<fik,ipi) / Wfe(£M(£>*) d £- 

J 



7 



Thus, a spectral decomposition of operator in (LRD) yields the following infinite-dimensional system 
on IV; of decoupled nth order equations: 



%k — A k x k , k — 0, 1, ... , 



(10) 



with x k (t) := [ x ljk (t) 



A h := 



b 2 - ttafc 







b n -i - a n i 



and ctift := Oj — Cji/fc = Oj + Ci(kir) 2 > 0. Based on [9j[T5] we conclude that each A k is Hurwitz if ([2]) 
holds. Therefore, each subsystem in fj 101) is exponentially stable and there exist Pk = P k > such that 

A\P k + P k A k = -I, k = 0,1, 

Now, since A is the infinitesimal generator of the following C -semigroup: 

oo 

T(tM0) := T(t)^,0) = ^ e A ^x k (0)ip k (0, 



k = 



we have 



) rao 00 / poo \ 

HT(^(o)n 2 dt : = y o (T(t)m,T(t)m) & = J2 x k(°)[j o e^v^j^o) 

oo 

= 5>£(o)p fc x fc (o). 



fc = 



We will show the exponential stability of the C -semigroup Tit) on -L^O, 1) by establishing convergence 
of the infinite sum X)ST=o x fc W^^CO) l° r eacn { x fc(0)}fceN € 1% [HI Lemma 5.1.2]. Let s m denote the 
mth partial sum, i.e. 



s m := J2xl(0)P k x k (0). 



k = 

For I < m we have 

m m 

\s m ~ Sl\ = 4(0)PkX k (0) < II^IHI^WII 2 - 

k=l+l k=l+l 

Now, we represent A k , for k ^ 0, as 

A k = k 2 (F + (l/k 2 )A ), F := -vr 2 diag{[ c x ■■■ c n ]} < 
and use perturbation analysis to express P k as 



(11) 



(12) 



k 2 



Pk = 77 + 79*i + 7^2 + 



A; 2 



fc 4 



1 1 

k 2 ^ k 2 i 3 

j = 



8 



where 

F V + V F = -I, F Vj + VjF = -(^-i + V*-i4>) (13) 
with j S N. Solution to (fl~3|) is determined by 



/•oo 

V = -{1/2)Fq\ Vj = / ^{AlV^ + V^A^dt 

Jo 

which can be used to obtain 

\\V \\ = l/(27T 2 C min ) 

\\VjW < ||F ||(2Po||||Vo||) J , j G N 

NT/ II 00 

Clearly, for k 2 > 2 ||-Ao|| ||^o|| the geometric series in the last inequality converges. This immediately gives 
the following upper bound for ||Pfe||: 

141 < m 



e - 2\\a,\\ II v&ir 

and inequality in (|12p simplifies to 

I I < ll^oll v 11 

|Sm - (/ + l) 2 -2p ||||Vo|| fc ^ +1 l|Xfc 

Hence, for each {£fc(0)}fc e N o G 2£ partial sum (fTTj) represents a Cauchy sequence which guarantees conver- 
gence of 

Sfc°=o and consequently 



DC 



||r(t)^(o)|| 2 dt < 00, w(o) g p(^). 







Since T>(A) is dense in ^(0, 1), by an argument as in 061 p. 51] this inequality can be extended to all 
ip(0) G L^iO, 1) which implies exponential stability of T(t) [471 Lemma 5.1.2]. 



Theorem 1 The C - semigroup T(t) generated by operator (LRD) (DM) on L^O, 1) is exponentially stable 
if the secant criterion (T2P is satisfied. 



3.2 The existence of a decoupled quadratic Lyapunov function 

The following theorem extends the diagonal stability result of [30] to PDEs with r spatial coordinates: 



Theorem 2 For system (LRD) (DM) there exist a decoupled quadratic Lyapunov function 



V{$) := (V,£>V>) = Y.diiiPi,^}, di > 0, (14) 



i = i 



that establishes exponential stability on LQiSl) if and only if holds. 
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Proof. We prove the theorem for a system given by 

ipt = M ■= CAip + A ^, 
where C := diag{[ c\ ••• c n ]} > 0, and 



(15) 



- 1 





... o 


-7i 


72 


- 1 










73 


-1 ■■. 










7n 


- 1 



(16) 



This is because all operators of the form (LRD) can be obtained by acting on Aq from the left with a 
diagonal matrix which does not change the existence of a decoupled quadratic Lyapunov function. We 
will prove that the secant criterion |(C3)| is both necessary and sufficient for the existence of a decoupled 
quadratic Lyapunov function. 

Necessity: Suppose that there exist a Lyapunov function of the form (|14|) that establishes exponential 
stability of (|15p . The derivative of (|14p along the solutions of (|15p is given by 



dV(i/>) 
dt 



{ipt,Dil>) + (^,£tyt) = (CAtp + Ao^D^) + (iP,DCAiP + DA i/;) 

n 

-2^c^(V^,VVi) + {i/>,(A^D + DAo)i>) < (ip, (AqD + DA )ij)' 



i = i 



where we have used Green's integral identity [19] with ip satisfying the Neumann boundary conditions on 
<9f2, and the fact that C and D commute. The exponential stability of (|15p and the above expression for 
dV(ip)/dt imply that Aq is Hurwitz. But |(C3)"| is a necessary condition for a matrix Aq with equal diagonal 
entries to be Hurwitz [9]. 

Sufficiency: Suppose that |(C3)| holds. Following [30] we define: 



r := (71 •••In) 1 > 0, T := diagjl, — , 

and differentiate (|14p along the solutions of (|15p to obtain 

dV(ip) 



i+l 72 



• • • 7n ) 

.n-1 / 



D := r - 



dt 



(I'M) 



If [(C3)1 holds then Q = Q T is a positive definite matrix [30] 

q := -(^o£> + da ) = -r-^rl^r- 1 + r- 1 A r)r- 1 > o 

and hence 



dV{jj) 
dt 



< -At, 



where A m i n ((5) > denotes the smallest eigenvalue of Q. Upon integration, we get 

< (m,Dm) < (V(0),^(0)) - A min (Q) TllrWWfdT 
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which yields 



^minlVJ 



(V»(o),dv(o)) , v< > o, v^(o) g p(^). 



Since £>(^4) is dense in L^O), the last inequality can be extended to all ip(0) G LV;(p) [361III]- Thus, for 
every ip(0) G Lg(fi) there is ^ := (^(0),D^(0)} /A mm (Q) > such that 



||T(i)V(0)|| 2 dr < N „ 



which proves the exponential stability of T(t) [471 Lemma 5.1.2]. 



Remark 1 The exponential stability ofT(t) in Theorem[J\ can be also established using a Lyapunov based 
approach with 

V(il>) = (i/>,Di/>), D := r- 2 diag{[ l/oi ••• l/a n ]}. 
However, the proof of Theorem [7] is of independent interest because of the explicit construction of the 



C -semigroup and block- diagonalization of operator (LRD) (DM) 



4 Extension to nonlinear reaction-diffusion equations 



We next show global asymptotic stability of the origin of the nonlinear distributed system (RD)tf(DM) 
This result holds in the LV^iQ) sense under the following assumption: 



Assumption 1 The functions fi(-), gi(-), and hi(-) in |(RD)| are continuously differentiable. Moreover, 
the functions fi(-) and gi{-) satisfy (CI) (C3), the functions hi(-) are positive, and the functions gi(-) are 
nondecreasing, i.e. 

K > 0, g ia := dgi/da > 0, Va G K. (C5) 

A new ingredient in Assumption [1] compared to the properties of /j(-) and gi{-) in ([5]) is a nondecreasing 
assumption on the functions gi(-). This additional assumption provides convexity of the Lyapunov function, 
which is essential for establishing stability in the presence of linear diffusion terms. For nonlinear diffusion 
terms we also assume that each /ij(-) is a positive function. 



Theorem 3 Suppose that system (RD) (DM) satisfies Assumption [7J Consider the Lyapunov function 
candidate 



where the di's are defined as in Section^ and suppose that there exists some function a(-) of class /Coo 
such that 

Vty) > a(U\\), W G (C6) 



Then if) = is a globally asymptotically stable equilibrium point of (RD) (DM), in the L^i^l) sense. 
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Remark 2 (Well-posedness) Standard arguments (see, for example, [36,50,51]) can be used to establish 
that (RD) (DM) has a unique solution on [0, i max ) • The existence of a unique solution on the time interval 



[0, oo) follows from the asymptotic stability of the origin of (RD) (DM) 



Proof. We represent the ith subsystem of (RD)f(DM) by: 

At = V • (hi(ipi)V7pi) ~ fibPi) + Ui 
Vi = 9i(ipi) 

k Au = 



H 



on dtt. 



The derivative of 



</>«(£) 



gi(a)da d£ 



(17) 



Vi(tpi) '■= li 
along the solutions of Hi is determined by 

7i(9ihl>i),ii>it) = H (9i(ipi), V ■ (hi(lj)i) Vipi) - fi(lpi) + Ui) . 

in combination with the Neumann boundary conditions on ipi, can be used 



Vi 



Green's integral identity 
to obtain 



Vi 



li {9i^ipi, hiVtpi) ~ li (9i, fi) + li (9i,Ui) . 



Now, from (C5) we have hiQi a > 0. Using this property and the fact that —lifi(cr)gi(a) < —gf(a) 
(cf. |(C1)||(C2)D we arrive at 

V < -(gi,gi) + n(gi,ui} = - {yuyi) + n(yi,Ui). 

This upper bound on Vi and the following Lyapunov function candidate: 



Vty) := X>VKV>0 



i = l 



yield 



V < (yMlD + DA )y) < -\ m UQ)h\ 



-A, 



xW)Ei 



9i 



(18) 



i = i 



Since the dj's are defined as in Section[3J we have used the fact that Q = Q T := — D + DAq) represents 
a positive definite matrix (see the proof of Theorem [2]). 

Now, since V(tp) > a(||?/>||) for each tp £ L^fi), with a(-) 6 /Coo, for any e > there exist 5 > 
such that ||?/>(0)|| < 5 implies ||^(f)|| < e for all t > 0. This follows from positive invariance of the set 
^fc := G ^2 (^)> ^(VO < &}> & > 0, and continuity of Lyapunov function 1/ [33]. Furthermore, V(V') is 
a nonincreasing function of time bounded below by zero and, thus, there exists a limit of V(ijj(t)) as time 
goes to infinity. If this limit is positive then (CI) , (C6) , and (I18|) imply the existence of m > such that 
su Pt>o V(tp(t)) < —m. But then V(ip(t)) < V(ip(0)) — mt and V(ip(t)) will eventually become negative 
which contradicts nonnegativity of V(ijj(t)), for all t > 0. Therefore, both V(ip(t)) and converge 



asymptotically to zero. ^,From the radial unboundedness of V(ip) (cf. (C6) ) and the above analysis we 
conclude global asymptotic stability of the origin, in the L^iP) sense. I 
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Remark 3 The condition |(C6)| on V(ip) can be weakened by working on L"(J7), in which case Jensen's 
inequality, applied to ( \17\) . provides the desired estimate (see Appendix [Tip , This relaxation allows for 
inclusion of many relevant nonlinearities arising in biological applications; one such example is provided 
in Section® Using a similar argument to the one presented in Theorem^ the global asymptotic stability 
of the origin in the L'i(£l) sense can be established (with keeping in mind that, in this case, (u,v) denotes 
a symbol for u T (£) v (£) d£j. 



5 Stability analysis for a compartmental model 



An alternative to the partial differential equation representation |(RD)| is a compartmental model which 
divides the reaction into compartments that are individually homogeneous and well-mixed, and represents 
them with ordinary differential equations. Compartmental models are preferable in situations where reac- 
tions are separated by physical barriers such as cell and intracellular membranes which allow limited flow 
between the compartments [52]. Instead of the lumped model © we now consider m compartments where 
the dynamics of the j'th compartment, j = 2, • • • ,m — 1, are given by 



x 3,l 
Xj,n 



l^j-l,2{Xj-l,2 



Xj,l) 
Xj, 2 ) 



Xj+l,l) 
Xj+1,2) 



- 9n{xj,n) 
f2(Xj,2) +9l(xj,l) 



(CM) 



3,nj 



The functions fJ>ji(-) i = 1, 
and possess the property 



,n,j = 1, 



m 



f^j,n{Xj t ri fn{,Xj,n)~\~9n—l{xj,n~l)- 

- 1, represent the diffusion terms between the compartments 



For the first and last compartments j 



o-[ij,i{o) > °> Va G R. (C7) 
= 1 and j = m, respectively the first and the second terms in the 



right-hand side of (CM) must be dropped because xo,i and x m +i,i are not defined. 



In the absence of the diffusion terms, the dynamics of the compartments in (CM) are decoupled, and 



coincide with (|5|) which is shown in [30] to be globally asymptotically stable under the conditions |(Cl)f 
|(C4)[ The following theorem makes an additional assumption that the function gi{-) be nondecreasing and 
proves that global asymptotic stability is preserved in the presence of diffusion terms: 



Theorem 4 Consider the compartmental model (CM) , j = 1, . . . , m, where for j = 1 and j = m, respec 



tively the first and the second terms in the right-hand side of (CM) are to be interpreted as zero. If the 
functions /j(-) and gi(-) satisfy the conditions (CI) (C4) and if, further, gi(-) is a nondecreasing function 
and //j,i(-) is as in |(C7)| then the origin Xj t i = is globally asymptotically stable. 

Proof. We first introduce the notation 

Xj := [ Xj : i 
= [ 



x 



Xj+l, 



.1.11 



i T 



J 
J 



m 



x 
1. 



[xj 



(19) 



In the absence of the diffusion terms in (CM) , the reference [30] constructs a Lyapunov function of the 
form 

n. 

L 3,i 



V(x 



i=l 



dm 



9i(o-)da 



(20) 
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where di, i = 1, • • • , n, are the diagonal entries of a matrix D obtained from ([6]) with A selected as in (1161) , 
and proves that its time derivative satisfies the estimate 



V(xj) < -e\\(gi(x jA ), - ■ ■ ,g n (xj, n ) 



for some e > 0. In the presence of the diffusion terms in (CM) , the estimate (|2ip becomes: 



V(xj) < -e\\{ gi {x hl ),--- ,g n {x hn ))\\ 2 + 



dV{ Xji 



while for j = 1: 
and for j = m: 



j = 2, • • ■ , m - 1 



V(xi) < -e||(ffi(a;ii),--- ,g n (xi, 



dV(xi) 
dx\ 



d Xj 



fn(xi - x 2 ) 



(21) 



(22) 



(23) 



V{x m ) < -e\\(gi(x m {),■■■ , g n {xm,n))\\ 2 + ' 



8x r 



Then the Lyapunov function 



v(x) = *£v( Xj ) 



satisfies 



V{x) < -e 



Substituting ([19]) and 



8Xn 



dl7iffi( 



(24) 



jj - Xj+i). (25) 



(26) 



which is obtained from f)20f) . we get 

m— 1 



<9.r 



'i+i 



m— 1 n 



j=i i=i 



Because gi(-) is a nondecreasing function by assumption, we note that \gi(xj t i) — <7i(:Cj+i,i)] possesses the 
same sign as (a;^ — Xj+ij). We next recall from property |(C7)| that fij^xj^ — also possesses the 

same sign as (xj^ — and, thus, 



\j9i( x j,i) ~ 9i(Xj+l,i)] VjA x j,i - x j+l,i) > 

which, according to (f2"7|) and ([25]) . implies 



V(x)<-e^||( 9l (x,. 



9n (Xj t r 



(28) 



(29) 



Because the Lyapunov function V(x) is proper from property | ( C4) | and because the right-hand side of (|29p 
is negative definite from property (CI) , we conclude that the origin x = is globally asymptotically stable. 
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Remark 4 Theorems^ and\^ both rely on the assumption that <?j(-) is nondecreasing, which translates to 
the convexity of the Lyapunov functions ( \17\) and ( \20\) . A similar convex Lyapunov function assumption 
has been employed in to preserve asymptotic stability in the presence of diffusion terms. Unlike the local 
result in JJ8j , however, in this paper we have established global asymptotic stability and allowed nonlinear 
diffusion terms by exploiting the specific structure of the system. 

6 An example 

We illustrate our main results with the analysis of a negative feedback loop around a simple MAPK cas- 
cade model. As described in the introduction, MAPK cascades are functional modules, highly conserved 
throughout evolution and across species, which mediate the transmission of signals generated by receptor 
activation into diverse biochemical and physiological responses involving cell cycle regulation, gene ex- 
pression, cellular metabolism, stress responses, and other functions. The control of MAPK and similar 
kinase cascades by therapeutic intervention is being investigated as a target for drugs, particularly in the 
areas of cancer and inflammation [53] . Several MAPK cascades have been found in yeast |54j and at least 
a dozen in mammalian cells [55j . and much effort is directed to the understanding of their dynamical 
behavior [28l[56l[57] . 

There are many models of MAPK cascades, with varying complexity. The simplest class of models |58| 
I18j . using quasi-steady state approximations for enzymatic mechanisms and a single phosphorylation site, 
involves a chain of three subsystems: 

hxi | u di(l - xi) 
ci+xi ei + (1 — x\) 
b 2 x 2 d 2 {l-x 2 ) 

-4- j;^ . 

c 2 + x 2 e 2 + (1 - x 2 ) 

b 3 x 3 d 3 (l-x 3 ) 

— T + x 2 —7z r 

C3 + X3 e 3 + (l-x 3 ) 

where u is an input and x 3 is seen as an output. The variables X{ denote the "active" forms of each of 
three proteins, and the terms 1 — x.; indicate the inactive forms of the respective proteins (after nondimen- 
sionalizing and assuming that the total concentration of each of the proteins = 1). For example, the term 
#1^2(1 — x 2 )/{e 2 + (1 — x 2 )) indicates the rate at which the inactive form of the second protein is being 
converted to active form. This rate is proportional to the concentration of the active form of the protein 
x\, which facilitates the conversion. Similarly, active x 2 facilitates the activation of the third protein. The 
first term in each of the right-hand sides models the inactivation of the respective protein, a mechanism 
that proceeds at a rate that is independent of the activation process. The saturated form of the nonlin- 
earities reflects the assumption that reactions are rate limited by resources such as the amount of enzymes 
available (an assumption that is not always valid). For this model, Kholodenko proposed in [18] the study 
of inhibitory feedback from the last to the first element, mathematically represented by a feedback law 
u = fi/(l + kx 3 ). See [18] for a description of the physical mechanism (an inhibitory phosphorylation of 
"SOS" protein, upstream of the system, by the last protein, p42/p44 MAPK or ERK) that might produce 
this inhibition. 

Linearizing the system about an equilibrium, there results a linear system to which one may apply the 
secant condition [18]. A linear model also arises when considering weakly activated pathways, the behavior 



xx = - 
x 2 = - 
x 3 = - 
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of the pathway when there is only a low level of kinase phosphorylation. In this case, one assumes that the 
inactive forms dominate: 1 — x\ ~ 1; this is the analysis in [58] and [59]. An intermediate case would be 
that in which activations are weak but the coefficients are small enough that the negative terms in the 
above equations cannot be replaced by linear functions; in that case, we are lead to equations as follows, 
for the closed-loop system: 

b\X\ fj, 



x\ = h 



ci + x\ 1 + kx 



3 

x 2 = h d 2 x\ (30) 

c 2 + x 2 

b 3 x 3 , 

X 3 = ; h d 3 X 2 

C3 + £3 

Both linearizations, as well as this nonlinear system, can be analyzed using our techniques. For our 
simulations we pick the nonlinear model, as it is more interesting. Denoting by x the equilibrium of (|30p 
and introducing the shifted variable x = x — x, we represent system (|30p as in ((5]) with 

fi(%i) = ~ ~~ _ j i = 1] 2, 3, 
C% ~r Xi C% "t X{ 

9i{xi) = d i+1 Xi, i = l,2, 



93(x 3 ) 



1 + kx 3 1 + kx 3 



We first note that condition (CI) is satisfied because fi(xi) and gi(xi), i = 1,2,3, are strictly increasing 
functions. Next, we recall from [30^ Section 6] that a set of gains 7^, i = 1,2,3, that do not depend on 
the specific location of x can be obtained by evaluating the maximum value of the slope ratio g\j f[ in the 
interval X{ £ [0, 1] in which Xi evolves. Upon trivial calculations we obtain: 

d i+1 { Ci + l) 2 kfi 2 (c 3 + l) 2 

li = 7 > *= 1,2 73 = 7 — max{c 3 , }. 

bid 63C3 (1 + k) 1 

We pick the parameters 6j = Cj = 1, i = 1, 2, 3, k = 1, d 2 = d 3 = fx = 0.4, which satisfy the secant 
condition (C3). (Although we chose the parameters to be 0(1) as in [58], we do not claim that these are 
physiologically realistic, nor are the diffusion constants that we pick below. We are merely interested in 
illustrating the theoretical results.) Adding diffusion terms with coefficients h\ = h 2 = h 3 = 0.001, we 
obtain the reaction-diffusion equations 



ht = 0.001 - — ^— + 



0.4 



1 + 01 1 + 03 

b 2t = 0.001 <j) 2 ^ - —p-r + °- 4< ^i ( EX ) 



1 + 



>2 



<j) 3t = 0.001 35C - + O.40 2 

1 + 03 

with the Neumann boundary conditions, 0^(0, t) = 0^(1, t) = 0. System [(EX)| can be brought to the 
form |(RD)1 using the following coordinate transformation: -0 := 0—0, where = [ 0.5501 0.2821 0.1272 
denotes the equilibrium point of |(EX)[ Asymptotic convergence of 0(x,t) to is illustrated in Fig. [3l A 
spatial discretization of the diffusion operator with Neumann boundary conditions is obtained using a 
Matlab Differentiation Matrix Suite |60|. 



T 
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Y>i(e,*) 



MM) 




Figure 3: Plots of := &(£,t) - 0, i = 1,2,3, for system |(EX)| with <£(f,0) = 

[ 16£ 2 (l-£ 2 ) 2 5 + costtC 2] T . 

7 Concluding remarks 

We identify a class of systems with a cyclic interconnection structure in which addition of diffusion does 
not have a destabilizing effect. For these systems, we demonstrate global stability if the "secant" criterion 
is satisfied. In the linear case, we show that the secant condition is necessary and sufficient for the existence 
of a decoupled Lyapunov function, which extends the diagonal stability result [30] to spatially distributed 
systems. For reaction-diffusion equations with nondecreasing coupling nonlinearities, we establish global 
asymptotic stability of the origin. Under some fairly mild assumptions, we also allow for nonlinear diffusion 
terms by exploiting the specific structure of the system. 



.1 Relaxation of condition (C6) 

Let us represent V^ipi) in (fT7|) by 



Vi(i/Ji) := 7i / Pi(ipi(0) d£> Pi( s ) : = / 9i(cr)da. 
Jq Jo 

and let O p (respectively, f2 m ) denote the set of points in Q where is positive (respectively, negative), 
i.e. 

n p -.= u e n, MO > 0} , n rn : = {£ e n, ^(£) < 0} . 

Then, Vi(ijji) can be rewritten as 



where 



- - /. 



Pi P (s) := / ffi(cr) da, p im {s) := I ft(er) da, s > 0. 
Jo Jo 

We observe that the first two derivatives of the functions pi p and pi m , respectively, satisfy 

{P' ip {s) = 9i(s) > 0, p'l p {s) = fas) > 0} 
{PiM = -9i(s) > 0, pUs) = g' t {-s) > 0} 
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which implies that both these functions are of class /Coo and convex. Using convexity, we may apply 
Jensen's inequality [39] to obtain 

Vi(lpi) > 7i (m>|pi p (||V>j||i/m>|) + l^mbimCIIV'illl/l^ml)) 

where \Q r \ denotes the measure of set r , and \\ipi\\i is the L\(0, l)-norm of fy. Since pi p and pi m are 
/Coo-functions we conclude that condition | ( C6) | on V{ip) always holds if the underlying state-space is L\ (U) 
(that is, there exists some function a(-) of class /Coo such that V(ip) > a(||^||i), Vf/) G L™(f2)). 
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